VMIproc class demo (image processing)

25/11/20 PH

Build on VMI class for VMI data analysis:

Implemented:

  • Plotting
  • Inversion
  • Uncertainties (bootstrap protocol)

To do:

  • Masking
  • Covariance

Background:

In [1]:
# # Memory profiler - OPTIONAL for testing
# # https://timothymonteath.com/articles/monitoring_memory_usage/
# %load_ext memory_profiler
# %memit


# # Quick hack for slow internet connection!
# %autosave 36000
In [2]:
# Standard imports
from pathlib import Path
import sys

import numpy as np
import xarray as xr

# Import class from file
sys.path.append(r'/cds/home/p/phockett/dev/tmo-dev/tmo')
import tmoDataBase as tb
# from vmi import VMI  # VMI class
from inversion import VMIproc  # VMI processing class

tb.setPlotDefaults()  # Set plot defaults (optional)

Load data

Set dataset

In [3]:
# Set main data path
baseDir = Path('/reg/data/ana15/tmo/tmolw0618/scratch/preproc/v2')

# Setup class object and which runs to read.
data = VMIproc(fileBase = baseDir, runList = range(89,97+1), fileSchema='_preproc_elecv2')

# The class has a bunch of dictionaries holding info on the data
# data.runs['files']

# Read in the files with h5py
# There is very basic checking implemented currently, just to see if 'energies' is present.
data.readFiles()
*** WARNING: key 90 missing energies data, will be skipped.
*** WARNING: key 92 missing energies data, will be skipped.
*** WARNING: key 94 missing energies data, will be skipped.
*** WARNING: key 95 missing energies data, will be skipped.
*** WARNING: key 96 missing energies data, will be skipped.
Read 9 files.
Good datasets: [89, 91, 93, 97]
Invalid datasets: [90, 92, 94, 95, 96]

Generate VMI images

Note this defaults to electron coords ('xc','yc').

In [4]:
data.genVMIXmulti()

# Check an image with Xarray plotter (Matplotlib)
data.imgStack['signal'].sel(run=89).pipe(np.log10).plot.imshow()
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/cds/home/p/phockett/.conda/envs/ps-4.0.8-local-clone/lib/python3.7/site-packages/xarray/core/computation.py:700: RuntimeWarning: divide by zero encountered in log10
  result_data = func(*input_data)
Out[4]:
<matplotlib.image.AxesImage at 0x7f46d741db90>

Set centre point

This is set as native pixel value, or assumed to be the max point in the image.

In [5]:
data.setCentre([553, 546])  # Set explictly
# data.setCentre()  # Default case - set as max ind.
data.checkCentre()
WARNING:param.OverlayPlot04820: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
In [6]:
# All the coord parameters are set in `centreInds`
# Note here that `cList` gives the array indicies for the centre point, based on the pixel coords set.
data.centreInds
Out[6]:
{'input': [553, 546],
 'dims': ['yc', 'xc'],
 'yc': {'coord': array(553.5), 'index': array([554])},
 'xc': {'coord': array(546.5), 'index': array([547])},
 'cList': [554, 547],
 'dMax': [493.0, 500.0]}

Basic processing

See the previous demo for info, here we'll just generate subtracted images.

In [7]:
# Subtraction, and set as new image set
data.imgStack['sub'] = data.imgStack['signal'] - data.imgStack['bg']
data.imgStack
Out[7]:
<xarray.Dataset>
Dimensions:  (run: 4, xc: 1048, yc: 1048)
Coordinates:
  * xc       (xc) float64 -0.5 0.5 1.5 2.5 ... 1.044e+03 1.046e+03 1.046e+03
  * yc       (yc) float64 -0.5 0.5 1.5 2.5 ... 1.044e+03 1.046e+03 1.046e+03
  * run      (run) int64 89 91 93 97
    norm     (run) int64 72318 73364 71472 36198
Data variables:
    signal   (xc, yc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    bg       (xc, yc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
    sub      (xc, yc, run) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
In [8]:
# View subtracted image stack
data.showImgSet(name='sub')
WARNING:param.RasterPlot07923: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.

Image inversion

Multiple methods to be supported... so far only Elio's cpBasex is implemented, and this needs a bit of setup.

In [9]:
# For cpbasex, set the basis path
basisPath = r'/reg/d/psdm/tmo/tmolw0618/scratch/results/tmo_ana/calc/G_r512_k128_l4.h5'

# Set also local version - pip installed version not working with loadG in testing (version mismatch issue, or something in import routine?)
pbasexPath = r'/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex' #/pbasex'

# Import module + basis functions
data.setMethod(pbasexPath = pbasexPath, basisPath = basisPath)
cpBasex import OK, checking for basis functions...
Found basis at /reg/d/psdm/tmo/tmolw0618/scratch/results/tmo_ana/calc/G_r512_k128_l4.h5.

Invert a set of images by setting filterSet (if not set, all images in self.imgStack will be processed). Optionally, set normalisation parameters for the images, the default is norm={'type':'max', 'scope':'global'}, which normalises to the global maximum (i.e. over all images).

Currently the options are:

type:

  • max
  • sum
  • raw

scope:

  • global
  • local
In [10]:
# Invert a dataset
data.inv(filterSet='signal')
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
In [11]:
# Full results are stored in a dictionay, self.proc[filterSet]
filterSet = 'signal'
data.proc[filterSet].keys()
Out[11]:
dict_keys(['fold', 'pbasex', 'xr'])
In [12]:
# ... where 'xr' contains the results in an Xarray
data.proc[filterSet]['xr']
Out[12]:
<xarray.DataArray 'IE' (BLM: 3, E: 512, run: 4)>
array([[[ 0.        ,  0.        ,  0.        ,  0.        ],
        [ 0.44566235,  0.43613563,  0.46423001,  0.43075811],
        [ 0.80048389,  0.77231224,  0.81993231,  0.79469168],
        ...,
        [ 0.05433236,  0.07269351,  0.08331476,  0.02396366],
        [ 0.05700265,  0.07613888,  0.08694888,  0.02535074],
        [ 0.05181977,  0.06914447,  0.07879847,  0.02314817]],

       [[        nan,         nan,         nan,         nan],
        [-2.03710948, -2.12816634, -2.17436883, -1.73430895],
        [-2.08275948, -2.20893407, -2.26650071, -1.74252256],
        ...,
        [ 0.68347204,  0.70466135,  0.7346497 ,  0.58669809],
        [ 0.84298988,  0.8666355 ,  0.91050975,  0.71909364],
        [ 0.92992677,  0.95506099,  1.00750427,  0.79006251]],

       [[        nan,         nan,         nan,         nan],
        [ 1.29333063,  1.41649754,  1.42861856,  1.22919167],
        [ 1.32871331,  1.4685956 ,  1.49439536,  1.22985456],
        ...,
        [-0.40262417, -0.44203047, -0.41306393, -0.39276484],
        [-0.56153601, -0.58910583, -0.58223553, -0.55475955],
        [-0.64865654, -0.66912341, -0.67599561, -0.64283519]]])
Coordinates:
  * BLM      (BLM) int64 0 2 4
  * E        (E) float64 0.0 0.000359 0.001436 0.003231 ... 93.01 93.38 93.74
  * run      (run) int64 89 91 93 97
    XS       (E, run) float64 0.0 0.0 0.0 0.0 30.82 ... 3.583 4.781 5.449 1.601
    mask     (E, run) bool False False False False False ... True True True True
    norm     (run) float64 59.72 60.73 69.15 62.71
Attributes:
    normType:  {'type': 'max', 'scope': 'global'}
In [13]:
# For rapid plotting per run, use self.plotSpectra()
data.plotSpectra()
In [14]:
# This can also be set to return a HoloMap for more flexibility
hmap = data.plotSpectra(returnMap = True)
hmap.overlay('BLM').layout('run').cols(1)
Out[14]:

Inversion details

Various other stages and results are also logged to the data strucure...

Images are set to 'filterSet-TYPE', where TYPE is:

  • 'fold': Quadrant filter result if quadrant filter is set. These are useful for checking centre point is OK. (TODO: implement quadrant.unfoldQuadrant for full symmetrized images.)
  • 'fit': Reconstructed/fitted images.
  • 'inv': Inverted images

These can be accessed and plotted with the usual functionality.

NOTE: there is currently no default radial masking set, so the default plot scaling will likely be off - adjust with the histogram.

In [15]:
# Quad filter 
data.showImg(name='signal-fold')
WARNING:param.RasterPlot18525: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
In [16]:
data.showImg(name='signal-inv')
WARNING:param.RasterPlot18875: Due to internal constraints, when aspect and width/height is set, the bokeh backend uses those values as frame_width/frame_height instead. This ensures the aspect is respected, but means that the plot might be slightly larger than anticipated. Set the frame_width/frame_height explicitly to suppress this warning.
In [17]:
# Raw cpbasex output data
data.proc['signal']['pbasex'].keys()
Out[17]:
dict_keys(['E', 'IE', 'betas', 'c', 'fit', 'inv'])

Uncertainties

A basis Poissonian bootstrapping routine is implemented for estimating uncertainties. This operates on the event data, and allows error propagation through the image generation & analysis, but is (consequently) a little slow. (For more on the method employed, see the wikipedeia page#Poisson_bootstrap).)

Note this routine is a little crude at the moment! It will run image generation for all filterSets, run image subtraction, then image inversion & stats for all filterSets, or just the listed element.

In [18]:
# Run bootstrapping & inversion routine, 5 iterations, and calculate stats for the subtracted data only.
data.bootstrapInv(N = 5, filterSet = 'sub')
*** Bootstrap run for N=5, lambda=1.0, filterSet=sub
Running set 1 of 5
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
Running set 2 of 5
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
Running set 3 of 5
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
Running set 4 of 5
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
Running set 5 of 5
Generating VMI images for filters: signal
Generating VMI images for filters: bg
/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex/pbasex.py:31: RuntimeWarning: invalid value encountered in true_divide
  betas = IEB[:,nim:].reshape(nx,nl-1,nim)/IE[:,None,:]
In [19]:
# Stats are output to a dataset, with variables per iteration
data.stats
Out[19]:
<xarray.Dataset>
Dimensions:  (BLM: 3, E: 512, run: 4)
Coordinates:
  * BLM      (BLM) int64 0 2 4
  * E        (E) float64 0.0 0.000359 0.001436 0.003231 ... 93.01 93.38 93.74
  * run      (run) int64 89 91 93 97
    XS       (E, run) float64 0.0 0.0 0.0 0.0 30.6 ... 3.249 4.362 4.946 1.439
    mask     (E, run) bool False False False False False ... True True True True
    norm     (run) float64 59.44 55.68 63.57 63.78
Data variables:
    0        (BLM, E, run) float64 0.0 0.0 0.0 0.0 ... -0.6575 -0.6655 -0.5502
    1        (BLM, E, run) float64 0.0 0.0 0.0 0.0 ... -0.6406 -0.6344 -0.5233
    2        (BLM, E, run) float64 0.0 0.0 0.0 0.0 ... -0.6264 -0.6143 -0.5827
    3        (BLM, E, run) float64 0.0 0.0 0.0 0.0 ... -0.6345 -0.6562 -0.6134
    4        (BLM, E, run) float64 0.0 0.0 0.0 0.0 ... -0.6961 -0.6839 -0.571
In [20]:
# Similarly to plotSpectra(), there is a plotUncertainties() function for this data
data.plotUncertainties()
In [21]:
# As before, an object can be returned for further plotting options
hmap = data.plotUncertainties(returnMap=True)
hmap.overlay('BLM') #.cols(1)
Out[21]:
In [22]:
hmap.overlay('run').layout('BLM').cols(1)
Out[22]:

Versions

In [23]:
import scooby
scooby.Report(additional=['holoviews', 'xarray'])
Out[23]:
Fri Nov 27 12:55:15 2020 PST
OS Linux CPU(s) 16 Machine x86_64
Architecture 64bit RAM 125.6 GB Environment Jupyter
Python 3.7.8 | packaged by conda-forge | (default, Jul 31 2020, 02:25:08) [GCC 7.5.0]
holoviews 1.13.5 xarray 0.16.1 numpy 1.19.4
scipy 1.5.3 IPython 7.19.0 matplotlib 3.3.2
scooby 0.5.6
In [24]:
# tmo-dev class version
data.__version__
Out[24]:
'0.0.1'